

* Open log
capture log close
log using "Data analysis\cbld-analysis01-fig04.log", replace text


* *************************************************
* Fig. 4. Observed versus predicted consensus rates
* *************************************************

* Programme:	cbld-analysis01-fig04.do
* Project:		Council coalition building
* Author:		Frank Haege, Department of Politics and Administration, University of Limerick
* Contact:		frank.haege@ul.ie

* Description
*************
* This do-file generates a figure comparing the observed consensus rates with the rates predicted
* by the model and the rates expected by chance. The figure also provides PRE statistics for the 
* predicted consensus rates for individual membership size periods and for the entire study period. 


* Set up Stata
version 11
clear all
macro drop _all
set linesize 80
set more off

* Load oberved consensus rates data
use "Experiments\Other data\cbld-management05a.dta", clear

* Rescale consensus rate variable
generate oconsensus = consensus/100

* Drop redundant variables
drop consensus

* Merge observed consensus rates data with predicted consensus rates data
merge m:m rulemembersize using "Experiments\Experiment02\cbld-management03-exp02a.dta"
drop _merge

* Drop simulation results for 27 member states (under Nice and Lisbon)
drop if year == .

* Calculate squared deviations
generate mdev = (oconsensus - consensus)^2
generate ndev = (oconsensus - econsensus)^2

* Generate sum of square
egen mss = sum(mdev), by(nostates)
egen tss = sum(ndev), by(nostates)
egen tmss = sum(mdev)
egen ttss = sum(ndev)
drop oconsensus mdev ndev

* Generate R-squared
generate r2 = (tss - mss) / tss
generate tr2 = (ttss - tmss) / ttss

* Collapse data set by number of member states
collapse (mean) r2=r2 tr2=tr2 consensus=consensus se=se lowerci=lowerci upperci=upperci /*
	*/ econsensus=econsensus ese=ese elowerci=elowerci eupperci=eupperci mss=mss tss=tss /*
	*/ tmss=tmss ttss=ttss, by(rulemembersize)

* Label variables
label var rulemembersize "Rule-membership period"
label var r2 "Prop. reduction of error for regime periods"
label var tr2 "Prop. reduction of error for entire period"
label var consensus "Expected consensus (model)"
label var se "Standard error (model)"
label var lowerci "CI lower bound (model)"
label var upperci "CI upper bound (model)"
label var econsensus "Expected consensus (null)"
label var ese "Standard error (null)"
label var elowerci "CI lower bound (null)"
label var eupperci "CI upper bound (null)"	
label var mss "Model sum of squares (within membership size)"
label var tss "Total sum of squares (within membership size)"
label var tmss "Model sum of squares (entire period)"
label var ttss "Total sum of squares (entire period)"

* Merge with observed consensus rates data
merge 1:m rulemembersize using "Experiments\Other data\cbld-management05b.dta"
drop _merge

* Generate variable for running number
generate no = _n
replace no = no - 3.2 if no > 3

* Merge with observed consensus rates data
merge m:m rulemembersize using "Experiments\Other data\cbld-management05c.dta"
drop _merge

* Update running number variable
replace no = _n if no == .
replace no = no - 6.1 if no > 6

* Label running number variable
label var no "Number of Member States"
label def nol 1 "12" 2 "15" 3 "25", replace
label val no nol

* Rescale variables
replace r2 = r2 *100
replace tr2 = tr2 *100
generate pconsensus = consensus *100
generate plowerci = lowerci * 100
generate pupperci = upperci * 100
generate peconsensus = econsensus * 100
generate pelowerci = elowerci * 100
generate peupperci = eupperci * 100

* Save R-squared measures in locals
sort rulemembersize
mkmat rulemembersize r2 tr2 mconsensus pconsensus peconsensus, matrix(r2)
matrix list r2
local r2_12 = string(round(r2[3,2], 1))
local r2_15 = string(round(r2[6,2], 1))
local r2_25 = string(round(r2[9,2], 1)) 
local tr2 = string(round(r2[3,3], 1))

* Plot model and null-model estimates against observed data
twoway (rcap pupperci plowerci no, /*
	*/ ylabel(0 20:100) xsize(6) /*
	*/ ytitle("Consensus Rate (Per Cent)") lwidth(medium) lcolor(grey) msize(medium) /*
	*/ xscale(range(0 3.5)) xlabel(1 2 3, val)) /*
	*/ (rcap peupperci pelowerci no, /*
	*/ lwidth(medium) lcolor(grey) msize(medium)) /*
	*/ (scatter peconsensus no, msymbol(s) mcolor(white)) /*
	*/ (scatter pconsensus no, msymbol(d) mcolor(white)) /*
	*/ (scatter mconsensus no, msymbol(dh) mcolor(white)) /*
	*/ (scatter consensus1994 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus1995 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus1996 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus1997 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus1998 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus1999 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus2000 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus2001 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus2002 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus2003 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus2004 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus2005 no, msymbol(oh) mcolor(white)) /*
	*/ (scatter consensus2006 no, msymbol(oh) mcolor(white)) /*
	*/ , legend(order(- "Predicted:" 4 "Computational Model" 3 "Null Model" /*
	*/ - " " - "Observed:" 5 "Mean" 6 "Individual Values" - " " - "Total {it:PRE}: `tr2'%") /*
	*/ cols(1) ring(1) position(3) region(lcolor(black))) scale(1.2) /*
	*/ text(10 .3 "{it:PRE} :") text(10 1 "`r2_12'%") text(10 2 "`r2_15'%") text(10 3 "`r2_25'%") /*
	*/ saving("Data analysis\Graphs\cbld-analysis01-fig01`i'.gph", replace)


log close
exit
